---
title: "Analysis of Medicare's Hospital Readmission Reduction Program"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: scroll
    source_code: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
library(RSocrata)
library(scales)
library(DT)
library(sp)
library(leaflet)
library(plotly)
library(tidyverse)

# Write a function to generate gauges for percentages
my_gauge <- function(x) {
  gauge(
    value = x, min = 0, max = 100,
    sectors = gaugeSectors(
      success = c(80, 100), warning = c(50, 79.99), danger = c(0, 49.99),
      colors = c("green", "yellow", "red")
    ),
    symbol = "%"
  )
}

# Read in hospital data
hosp_data <- read_rds("ghi_data.rds") %>%
  
  # Keep only the necessary variables
  select(
    provider_id, state, hospital_name, hospital_overall_rating, hospital_type,
    location.coordinates, ends_with("national_comparison")
  ) %>%
  
  # Clean up the variable names
  set_names(
    names(.) %>% str_replace_all("hospital_", "") %>% 
      str_replace_all("_national_comparison", "")
  ) %>%
  
  dplyr::rename(
    timeliness = "timeliness_of_care", 
    effectiveness = "effectiveness_of_care",
    imaging = "efficient_use_of_medical_imaging",
    experience = "patient_experience",
    safety = "safety_of_care"
  ) %>%
  
  # Replace all the values of "Not Available" with NA
  mutate_all(
    funs(replace(., . == "Not Available", NA_character_))
  ) %>%
  
  # Covert all categorical variables to factor
  mutate_at(
    vars(effectiveness:timeliness), 
    funs(
      str_to_title(.) %>%
      fct_relevel(
        "Above The National Average", "Same As The National Average",
        "Below The National Average"
      )
    )
  ) %>% 
  
  # Reorder the overall rating variable from 1 to 5 and create a status variable
  # indicating whether a hospital is monitored by the HRRP
  mutate(
    overall_rating = as.character(overall_rating) %>% 
      fct_relevel(as.character(1:5)),
    status = ifelse(
      str_detect(type, "^Acute") & !state == "MD", "HRRP", "Exempt")
  )

# Unlist the location data
hosp_data <- hosp_data %>%
  left_join(
    hosp_data %>% select(provider_id, location.coordinates) %>%
      drop_na(location.coordinates) %>%
      mutate(
        long = flatten_dbl(
          map(location.coordinates, ~ ifelse(!is.null(.x), .x[1], .x))
        ), 
        lat = flatten_dbl(
          map(location.coordinates, ~ ifelse(!is.null(.x), .x[2], .x))
        )
      ) %>%
      select(-location.coordinates),
    by = "provider_id"
  ) %>%
  select(-location.coordinates)


# Read in the program data
prog_data <- read_rds("hrrp_data.rds") %>%
  
  # Keep only necessary variables
  select(provider_id, measure_id, readm_ratio) %>%
  
  # Keep only the important parts of the measure ID label names and convert it
  # to a factor variable and convert the readmission ratio data to numeric
  mutate(
    measure_id = str_replace_all(measure_id, "^READM-30-", "") %>%
      str_replace_all("-HRRP", "") %>% as_factor,
    readm_ratio = as.numeric(readm_ratio)
  )

# Merge the data
hrrp_data <- hosp_data %>%
  left_join(
    prog_data %>%
      
      # Spread the column with the conditions across multiple columns using the
      # readmission ratio as the values
      spread(measure_id, readm_ratio) %>%
      
      # Create a variable indicating whether the given hospital has readmissions
      # rate data
      mutate(
        rprt = select(., -provider_id) %>% 
          apply(1, function(x) as.numeric(sum(!is.na(x)) > 0)),
        mn_rr = rowMeans(select(., -provider_id, -rprt), na.rm = TRUE)
      ),
    by = "provider_id"
  )
```

Main
================================================================================

Counts
--------------------------------------------------------------------------------

```{r gen_counts, echo=FALSE, include=FALSE}
# Generate sample counts, counts of hospitals reporting, total readmissions per
# hospital, and percent reporting
counts <- hrrp_data %>%
  group_by(status) %>%
  summarise(
    count = n(),
    rprt_count = sum(rprt, na.rm = T), 
    mean_rr = mean(mn_rr, na.rm = T)
  ) %>%
  bind_rows(
    hrrp_data %>%
      summarise(
        count = n(),
        rprt_count = sum(rprt, na.rm = T), 
        mean_rr = mean(mn_rr, na.rm = T)
      ) %>%
      mutate(status = "All")
  ) %>%
  ungroup %>%
  mutate(pct_rprt = scales::percent(rprt_count / count))
```

### Sample Size: All
```{r sample_all}
valueBox(counts %>% filter(status == "All") %>% pull(count))
```


### Percent Reorting: All
```{r pct_rprt_all}
my_gauge(counts %>% filter(status == "All") %>% pull(pct_rprt))
```

### Sample Size: HRRP
```{r sample_hrrp}
valueBox(counts %>% filter(status == "HRRP") %>% pull(count))
```


### Percent Reorting: HRRP
```{r pct_rprt_hrrp}
my_gauge(counts %>% filter(status == "HRRP") %>% pull(pct_rprt))
```

### Sample Size: Exempt
```{r sample_exempt}
valueBox(counts %>% filter(status == "Exempt") %>% pull(count))
```


### Percent Reorting: Exempt
```{r pct_rprt_exempt}
my_gauge(counts %>% filter(status == "Exempt") %>% pull(pct_rprt))
```


Violin Plots
--------------------------------------------------------------------------------

### Violin Plot: HRRP
```{r vio_plot_hrrp}
hosp_data %>%
  left_join(prog_data, by = "provider_id") %>%
  filter(status %in% "HRRP") %>%
  drop_na(measure_id, readm_ratio) %>%
  plot_ly(
    x = ~measure_id, y = ~readm_ratio, split = ~measure_id, type = 'violin',
    meanline = list(visible = TRUE)
  ) %>%
  layout(
    xaxis = list(title = "Condition"), 
    yaxis = list(title = "Readmission Ratio"),
    title = "Readmission Ratio Distributions Across Conditions: HRRP",
    showlegend = FALSE
  )
```

### Violin Plot: Exempt
```{r vio_plot_exempt}
hosp_data %>%
  left_join(prog_data, by = "provider_id") %>%
  filter(status %in% "Exempt") %>%
  drop_na(measure_id, readm_ratio) %>%
  plot_ly(
    x = ~measure_id, y = ~readm_ratio, split = ~measure_id, type = 'violin',
    meanline = list(visible = TRUE)
  ) %>%
  layout(
    xaxis = list(title = "Condition"), 
    yaxis = list(title = "Readmission Ratio"),
    title = "Readmission Ratio Distributions Across Conditions: Exempt",
    showlegend = FALSE
  )
```


Metric Comparison Plots
--------------------------------------------------------------------------------

```{r gen_comp_data, echo=FALSE, include=FALSE}
comp_plot_data <- hrrp_data %>%
  select(effectiveness:timeliness, Status = "status") %>%
  gather(Metric, nc, effectiveness:timeliness) %>%
  drop_na %>%
  mutate(
    Metric = str_to_title(Metric),
    nc = fct_relevel(
      nc, "Above The National Average", "Same As The National Average", 
      "Below The National Average"
    )
  ) %>%
  count(Status, Metric, nc) %>%
  group_by(Status, Metric) %>%
  mutate(Percentage = n * 100 / sum(n))
```

### Metric Comparison Plot: HRRP
```{r comp_plot_hrrp}
comp_plot_data %>% 
  filter(Status == "HRRP") %>% 
  plot_ly(x = ~Metric, y = ~Percentage, color = ~nc, type = 'bar') %>%
  layout(
    title = "Metric Comparisons of HRRP Hospitals",
    barmode = 'group',
    legend = list(orientation = 'h'),
    xaxis = list(title = "")
  )
```

### Metric Comparison Plot: Exempt
```{r comp_plot_exempt}
comp_plot_data %>% 
  filter(Status == "Exempt") %>% 
  plot_ly(x = ~Metric, y = ~Percentage, color = ~nc, type = 'bar') %>%
  layout(
    title = "Metric Comparisons of Exempt Hospitals",
    barmode = 'group',
    legend = list(orientation = 'h'),
    xaxis = list(title = "")
  )
```


Map
--------------------------------------------------------------------------------
```{r gen_geo_data, echo=FALSE, include=FALSE}
geo_data <- hrrp_data %>%
  filter(status %in% "HRRP") %>%
  select(long, lat, readm_ratio = "mn_rr", name) %>%
  drop_na()

cols <- colorQuantile(c("blue", "gray", "red"), geo_data$readm_ratio, n = 5)

geo_data <- SpatialPointsDataFrame(
  geo_data %>% select(long, lat), 
  geo_data %>% select(-c(long, lat))
)
```

### HRRP Map
```{r}
# Generate the map
leaflet(data = geo_data) %>%
  addTiles() %>%
  
  # Set the center of the map and the zoom level
  setView(-98.5795, 39.8283, zoom = 2.3) %>%
  
  # Add the markers colored based on whether the given hospital ratio's was high
  # or low and add a label that shows the hospital name and its readmission
  # ratio
  addCircleMarkers(
    radius = 1,
    color = ~cols(readm_ratio),
    opacity = 0.8,
    label = ~paste(name, "-", readm_ratio)
  ) %>%
  addLegend(
    "bottomright", pal = cols, values = ~readm_ratio,
    title = "Readmission Ratio Quintile", opacity = 1
  )
```


Data
================================================================================

### HRRP Data
```{r hrrp_table}
DT::datatable(hrrp_data, options = list(pageLength = 10, scrollX = TRUE))
```